from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
from ngsolve.bem import *

keys: Acoustic Scattering, Brakhage-Werner formulation, Helmholtz equation, Dirichlet data

Helmholtz Problem

We consider a Helmholtz problem with a plane wave as initial condition, i.e., \(m = e^{i\,\kappa \,x} \) and

\[\begin{split} \left\{ \begin{array}{rcl l} \Delta u + \kappa^2 \, u &=&0, \quad &\Omega \subset \mathbb R^3\,,\\ \gamma_0\, u &=& m, \quad &\Gamma = \partial \Omega\, .\end{array} \right. \end{split}\]

Combined field integral equations combine single and double layer integral operators, one simple option is the Brakhage-Werner formulation.

The solution is represented as

\[ u = i \,\kappa \,\mathrm{SL}(j) - \mathrm{DL}(j)\,, \]

where \(j\) solve the boundary integral equation $\( \left( \frac{1}{2} + \mathrm K + i \, \kappa \,\mathrm V \right) \, \mathrm j = \mathrm u_{in} \qquad \text{on} \, \Gamma = \partial \Omega \)$

The single layer and double layer operator are defined using the fundamental solution of the Helmholtz operator,

\[ G(x,y) = \frac{e^{i\kappa \|x-y\|}}{4\pi\,\|x-y\|}. \]

The layer potentials \( \mathrm SL\) and \(\mathrm DL\) have the same form as in the Laplace case, except that the Laplace kernel is replaced by the Helmholtz kernel \(G\).

kappa=10
order=4

screen = WorkPlane(Axes( (0,0,0), Z, X)).RectangleC(15,15).Face()

sp = Fuse(Sphere( (0,0,0), pi).faces)
screen.faces.name="screen"
sp.faces.name="sphere"
shape = Compound([screen,sp])

mesh = shape.GenerateMesh(maxh=5/kappa).Curve(order)
Draw (mesh);

fes_sphere = Compress(SurfaceL2(mesh, 
                                order=order, 
                                complex=True, 
                                definedon=mesh.Boundaries("sphere")))
u,v = fes_sphere.TnT()
fes_screen = Compress(SurfaceL2(mesh, order=order, 
                                dual_mapping=True, 
                                complex=True, 
                                definedon=mesh.Boundaries("screen")))
print ("ndof_sphere = ", fes_sphere.ndof, "ndof_screen =", fes_screen.ndof)

ndof_sphere =  17010 ndof_screen = 31110

with TaskManager():
    C = HelmholtzCF(u*ds("sphere"), kappa)*v*ds
    u,v  = fes_sphere.TnT()
    Id = BilinearForm(u*v*ds).Assemble()

lhs = 0.5 * Id.mat + C.mat
m = exp(1j * kappa * x) 
rhs = LinearForm(m*v*ds).Assemble()

gfu = GridFunction(fes_sphere)
pre = BilinearForm(u*v*ds, diagonal=True).Assemble().mat.Inverse()
with TaskManager():
    gfu.vec[:] = solvers.GMRes(A=lhs, b=rhs.vec, pre=pre, maxsteps=40, tol=1e-8)

GMRes iteration 1, residual = 60.52440579389444     

GMRes iteration 2, residual = 24.001710415059282     

GMRes iteration 3, residual = 10.81618663970256     

GMRes iteration 4, residual = 5.17850130118241     

GMRes iteration 5, residual = 2.549268448878951     

GMRes iteration 6, residual = 1.288805495688266     

GMRes iteration 7, residual = 0.6549832602491101     

GMRes iteration 8, residual = 0.3554164558330166     

GMRes iteration 9, residual = 0.19990866606606392     

GMRes iteration 10, residual = 0.10836869993844735     

GMRes iteration 11, residual = 0.059314453140396746     

GMRes iteration 12, residual = 0.032946081231833756     

GMRes iteration 13, residual = 0.013819139108246862     

GMRes iteration 14, residual = 0.007360642470274827     

GMRes iteration 15, residual = 0.0033138105281745843     

GMRes iteration 16, residual = 0.0019472886266694725     

GMRes iteration 17, residual = 0.001152960819267155     

GMRes iteration 18, residual = 0.0006865451125338632     

GMRes iteration 19, residual = 0.00040647668374400083     

GMRes iteration 20, residual = 0.00023525302450811517     

GMRes iteration 21, residual = 0.00012568335128557858     

GMRes iteration 22, residual = 7.27284483870596e-05     

GMRes iteration 23, residual = 4.282868528421723e-05     

GMRes iteration 24, residual = 2.238947095910078e-05     

GMRes iteration 25, residual = 1.1323654081077384e-05     

GMRes iteration 26, residual = 5.6817723197032995e-06     

GMRes iteration 27, residual = 2.815481215969824e-06     

GMRes iteration 28, residual = 1.3918984915180687e-06     

GMRes iteration 29, residual = 7.931632776571698e-07     

GMRes iteration 30, residual = 4.523061248529819e-07     

GMRes iteration 31, residual = 2.723734968721705e-07     

GMRes iteration 32, residual = 1.6096117668515718e-07     

GMRes iteration 33, residual = 9.476709863740349e-08     

GMRes iteration 34, residual = 5.214502615501129e-08     

GMRes iteration 35, residual = 2.9709338426800344e-08     

GMRes iteration 36, residual = 1.733994073439006e-08     

GMRes iteration 37, residual = 9.167410459121383e-09     

Draw (gfu, order=5, min=-1, max=1);

Postprocessing on screen

uscat = GridFunction(fes_screen)
with TaskManager():
    uscat.Set(1j*kappa*HelmholtzSL(u*ds("sphere"),kappa)(gfu) - HelmholtzDL(u*ds("sphere"), kappa)(gfu), \
              definedon=mesh.Boundaries("screen"))

print ("Scattered field")
Draw (uscat, mesh, min=-1,max=1, animate_complex=True, order=4);

Scattered field

uin = mesh.BoundaryCF( {"screen": m }, default=0)
print ("Total field")
Draw (uin+uscat, mesh, min=-1,max=1, animate_complex=True, order=4);

Total field